• Методы построения деревьев
• Методы кластеризации на основании расстояний
• Примеры для демонстрации и для заданий
• Кластерный анализ в R
• Качество кластеризации:
  • кофенетическая корреляция
  • ширина силуэта
  • поддержка ветвей
• Сопоставление деревьев: танглграммы

Вы сможете

• Выбирать подходящий метод агрегации (алгоритм кластеризации)
• Строить дендрограммы
• Оценивать качество кластеризации (Кофенетическая корреляция, ширина силуэта, поддержка ветвей)
• Сопоставлять дендрограммы, полученные разными способами, при помощи танглграмм

Методы кластеризации на основании расстояний (о них сегодня)

• Метод ближайшего соседа
• Метод отдаленного соседа
• Метод среднегруппового расстояния
• Метод Варда
• и т.д. и т.п.

Методы кластеризации на основании признаков

• Метод максимальной бережливости
• Метод максимального правдоподобия

И это еще далеко не все

Морфометрия черепов у волков в Скалистых горах и в Арктике (Jolicoeur, 1959)

Map from Jolicoeur 1959

Данные взяты из работы Morrison (1990):

• A — волки из Арктики (10 самцов, 6 самок)
• L — волки из Скалистых гор (6 самцов, 3 самки)

library(candisc)
data("Wolves")

dim(Wolves)

## [1] 25 12

colnames(Wolves)

##  [1] "group"    "location" "sex"      "x1"       "x2"       "x3"      
##  [7] "x4"       "x5"       "x6"       "x7"       "x8"       "x9"

head(rownames(Wolves))

## [1] "rmm1" "rmm2" "rmm3" "rmm4" "rmm5" "rmm6"

any(is.na(Wolves))

## [1] FALSE

table(Wolves$group)

## 
## ar:f ar:m rm:f rm:m 
##    6   10    3    6

data("iris")

dim(iris)

## [1] 150   5

colnames(iris)

## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 
## [5] "Species"

head(rownames(iris))

## [1] "1" "2" "3" "4" "5" "6"

# Делаем осмысленные имена строк
Species <- substr(iris$Species, 0, 2)
rownames(iris) <- make.unique(Species)
# Делаем случайную выборку для этой демонстрации
set.seed(191231)
ids <- sample(nrow(iris), 50)
siris <- iris[ids, ]

• Постройте ординацию nMDS данных о морфометрии волков и ирисов
• Оцените качество ординации
• Обоснуйте выбор коэффициента
• Раскрасьте точки на ординации волков в зависимости от географического происхождения (group), а на ординации ирисов — от вида (Species)

library(vegan, quietly = TRUE)
library(ggplot2); theme_set(theme_bw())
s_w <- scale(Wolves[, 4:ncol(Wolves)]) ## стандартизируем
ord_w <- metaMDS(comm = s_w, distance = "euclidean", autotransform = FALSE)

## Run 0 stress 0.101 
## Run 1 stress 0.143 
## Run 2 stress 0.101 
## ... Procrustes: rmse 0.0000114  max resid 0.0000343 
## ... Similar to previous best
## Run 3 stress 0.128 
## Run 4 stress 0.128 
## Run 5 stress 0.101 
## ... Procrustes: rmse 0.0000221  max resid 0.0000781 
## ... Similar to previous best
## Run 6 stress 0.138 
## Run 7 stress 0.124 
## Run 8 stress 0.101 
## ... Procrustes: rmse 0.000011  max resid 0.0000363 
## ... Similar to previous best
## Run 9 stress 0.127 
## Run 10 stress 0.147 
## Run 11 stress 0.155 
## Run 12 stress 0.135 
## Run 13 stress 0.137 
## Run 14 stress 0.127 
## Run 15 stress 0.155 
## Run 16 stress 0.101 
## ... Procrustes: rmse 0.00000238  max resid 0.0000052 
## ... Similar to previous best
## Run 17 stress 0.137 
## Run 18 stress 0.135 
## Run 19 stress 0.127 
## Run 20 stress 0.154 
## *** Solution reached

dfr_w <- data.frame(ord_w$points, Group = Wolves$group)
gg_w <- ggplot(dfr_w, aes(x = MDS1, y = MDS2)) + geom_point(aes(colour = Group))

ord_i <- metaMDS(comm = siris[, -5], distance = "euclidean", 
                 autotransform = FALSE)

## Run 0 stress 0.0272 
## Run 1 stress 0.0336 
## Run 2 stress 0.0272 
## ... Procrustes: rmse 0.000325  max resid 0.000901 
## ... Similar to previous best
## Run 3 stress 0.0319 
## Run 4 stress 0.034 
## Run 5 stress 0.0275 
## ... Procrustes: rmse 0.0125  max resid 0.0804 
## Run 6 stress 0.0362 
## Run 7 stress 0.0272 
## ... Procrustes: rmse 0.0000656  max resid 0.000199 
## ... Similar to previous best
## Run 8 stress 0.0314 
## Run 9 stress 0.0314 
## Run 10 stress 0.0321 
## Run 11 stress 0.0379 
## Run 12 stress 0.0358 
## Run 13 stress 0.0294 
## Run 14 stress 0.0319 
## Run 15 stress 0.0356 
## Run 16 stress 0.0362 
## Run 17 stress 0.0323 
## Run 18 stress 0.0272 
## ... Procrustes: rmse 0.0000111  max resid 0.0000599 
## ... Similar to previous best
## Run 19 stress 0.0289 
## Run 20 stress 0.0294 
## *** Solution reached

dfr_i <- data.frame(ord_i$points, Species = siris$Species)
gg_i <- ggplot(dfr_i, aes(x = MDS1, y = MDS2)) + geom_point(aes(colour = factor(Species)))

Результат кластеризации зависит от

• коэффициента сходства-различия
• от алгоритма кластеризации

= nearest neighbour = single linkage

Особенности:

• Может быть сложно интерпретировать, если нужны группы
  • объекты на дендрограмме часто не образуют четко разделенных групп
  • часто получаются цепочки кластеров (объекты присоединяются как бы по-одному)
• Хорош для выявления градиентов

= furthest neighbour = complete linkage

Особенности:

• На дендрограмме образуется много отдельных некрупных групп
• Хорош для поиска дискретных групп в данных

= UPGMA = Unweighted Pair Group Method with Arithmetic mean

Особенности:

• UPGMA и WUPGMС иногда могут приводить к инверсиям на дендрограммах

Инверсии на дендрограммах

= Ward’s Minimum Variance Clustering

Особенности:

• метод годится и для неевклидовых расстояний несмотря на то, что внутригрупповая дисперсия расстояний рассчитывается так, как будто это евклидовы расстояния

Давайте построим деревья при помощи нескольких алгоритмов кластеризации (по стандартизованным данным, с использованием Евклидова расстояния) и сравним их.

# Нам понадобится матрица расстояний
d <- dist(x = s_w, method = "euclidean")
# Пакеты для визуализации кластеризации
library(ape)
library(dendextend)

hc_single <- hclust(d, method = "single")

• И это все?
• Нет!

plot(hc_single)

ph_single <- as.phylo(hc_single)
plot(ph_single, type = "phylogram", cex = 0.7)
axisPhylo()

den_single <- as.dendrogram(hc_single)
plot(den_single)

hc_compl <- hclust(d, method = "complete")
ph_compl <- as.phylo(hc_compl)
plot(ph_compl, type = "phylogram", cex = 0.7)
axisPhylo()

den_compl <- as.dendrogram(hc_compl)
plot(den_compl)

hc_avg <- hclust(d, method = "average")
ph_avg <- as.phylo(hc_avg)
plot(ph_avg, type = "phylogram", cex = 0.7)
axisPhylo()

den_avg <- as.dendrogram(hc_avg)
plot(den_avg)

hc_w2 <-hclust(d, method = "ward.D2")
ph_w2 <- as.phylo(hc_w2)
plot(ph_w2, type = "phylogram", cex = 0.7)
axisPhylo()

den_w2 <- as.dendrogram(hc_w2)
plot(den_w2)

Кофенетическая корреляция - мера качества отображения многомерных данных на дендрограмме. Кофенетическую корреляцию можно рассчитать как Пирсоновскую корреляцию (обычную) между матрицами исходных и кофенетических расстояний между всеми парами объектов

Кофенетическое расстояние - расстояние между объектами на дендрограмме

Метод агрегации, который дает наибольшую кофенетическую корреляцию, дает кластеры, лучше всего отражающие исходные данные

Можно рассчитать при помощи функции из пакета ape

c_single <- cophenetic(ph_single)
c_compl <- cophenetic(ph_compl)
c_avg <- cophenetic(ph_avg)
c_w2 <- cophenetic(ph_w2)

cor(d, as.dist(c_single))

## [1] 0.565

cor(d, as.dist(c_compl))

## [1] 0.706

cor(d, as.dist(c_avg))

## [1] 0.745

cor(d, as.dist(c_w2))

## [1] 0.726

Оцените для данных об ирисах при помощи кофенетической корреляции качество кластеризаций, полученных разными методами.

Какой метод дает лучший результат?

• Можно выбрать число кластеров:
  • либо субъективно, на любом выбранном уровне (главное, чтобы кластеры были осмысленными и интерпретируемыми);
  • либо исходя из распределения расстояний ветвления.
• Можно оценить стабильность кластеризации при помощи бутстрепа.

Ширина силуэта \(s_i\) — мера степени принадлежности объекта \(i\) к кластеру. Это среднее расстояние от данного объекта до других объектов из того же кластера (\(\bar{d}_{i~within}\)), в сравнении со средним расстоянием до ближайшего кластера (\(\bar{d}_{i~to~nearest~cluster}\)).

\(s_i = \frac {\bar{d}_{i~to~nearest~cluster} - \bar{d}_{i~within}} {max\{\bar{d}_{i~to~nearest~cluster}, \bar{d}_{i~within}\}}\)

\(-1 \le s_i \le 1\) — чем больше \(s_i\), тем “лучше” объект принадлежит кластеру.

• Средняя ширина силуэта для всех объектов из кластера — оценивает, насколько “тесно” сгруппированы объекты.
• Средняя ширина силуэта по всем данным — оценивает общее качество классификации.
• Чем больше к 1, тем лучше. Если меньше 0.25, то можно сказать, что нет структуры.

Оценим ширину силуэта для 3 кластеров

library(cluster)
complete3 <- cutree(tree = hc_avg, k = 3)
plot(silhouette(x = complete3, dist = d), cex.names = 0.6)

“An approximately unbiased test of phylogenetic tree selection” (Shimodaria, 2002)

Этот тест использует специальный вариант бутстрепа — multiscale bootstrap. Мы не просто многократно берем бутстреп-выборки и оцениваем для них вероятность получения топологий (BP p-value), эти выборки еще и будут с разным числом объектов. По изменению BP при разных объемах выборки можно вычислить AU (approximately unbiased p-value).

library(pvclust)
# итераций должно быть 1000 и больше, здесь мало для скорости
cl_boot <- pvclust(t(s_w), method.hclust = "average", nboot = 100, 
                   method.dist = "euclidean", parallel = TRUE, iseed = 42)

## Creating a temporary cluster...done:
## socket cluster with 3 nodes on host 'localhost'
## Multiscale bootstrap... Done.

Обратите внимание на число итераций: nboot = 100 — это очень мало. На самом деле нужно 10000 или больше.

AU — approximately unbiased p-values (красный), BP — bootstrap p-values (зеленый)

plot(cl_boot)

# pvrect(cl_boot) # достоверные ветвления

График стандартных ошибок для AU p-value нужен, чтобы оценить точность оценки самих AU. Чем больше было бутстреп-итераций, тем точнее будет оценка AU.

seplot(cl_boot, cex = 0.5)

# print(cl_boot) # все значения

Два дерева (с непохожим ветвлением) выравнивают, вращая случайным образом ветви вокруг оснований. Итеративный алгоритм. Картина каждый раз разная.

set.seed(395)
untang_w <- untangle_step_rotate_2side(den_compl, den_w2, print_times = F)

# танглграмма
tanglegram(untang_w[[1]], untang_w[[2]],
           highlight_distinct_edges = FALSE,
           common_subtrees_color_lines = F,
           main = "Tanglegram",
           main_left = "Left tree",
           main_right = "Right tree",
           columns_width = c(8, 1, 8),
           margin_top = 3.2, margin_bottom = 2.5,
           margin_inner = 4, margin_outer = 0.5,
           lwd = 1.2, edge.lwd = 1.2, 
           lab.cex = 1, cex_main = 1)

Постройте танглграмму для данных о морфометрии ирисов из дендрограмм, полученных методом ближайшего соседа и методом Варда.

webpage <-"http://evolution.genetics.washington.edu/book/primates.dna"
primates.dna <- read.dna(webpage)
d_pri <- dist.dna(primates.dna)
hc_pri <- hclust(d_pri, method = "average")
ph_pri <- as.phylo(hc_pri)
plot(ph_pri)
axisPhylo()

• Результат кластеризации зависит не только от выбора коэффициента, но и от выбора алгоритма.
• Качество кластеризации можно оценить разными способами.
• Кластеризации, полученные разными методами, можно сравнить на танглграммах.

• Borcard, D., Gillet, F., Legendre, P., 2011. Numerical ecology with R. Springer.
• Legendre, P., Legendre, L., 2012. Numerical ecology. Elsevier.
• Quinn, G.G.P., Keough, M.J., 2002. Experimental design and data analysis for biologists. Cambridge University Press.

Как работает UPGMA можно посмотреть здесь:

• http://www.southampton.ac.uk/~re1u06/teaching/upgma/

Как считать поддержку ветвей (пакет + статья):

• pvclust: An R package for hierarchical clustering with p-values [WWW Document], n.d. URL http://www.sigmath.es.osaka-u.ac.jp/shimo-lab/prog/pvclust/ (accessed 11.7.14).

Для анализа молекулярных данных:

• Paradis, E., 2011. Analysis of Phylogenetics and Evolution with R. Springer.